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ABSTRACT 

A formulation is presented that describes the 
macrotheromechanical behavior of materials subjected 
to rapid thermal or mechanical loading such as occurs 
in most heat treatments. The equations are developed 
for Lagrangian, Eulerian, and intermediary kinematic 
descriptions and are independent of the constitutive 
laws and the equation of state; they can be solved 
numerically for a specified material and boundary 
conditions. The coupled transport effects between 
dissipation and energy are included. The conventional 
linearized stability approach indicates the numerical 
procedure to be stable, with certain restrictions on 
the time step size. 
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• time derivative 
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j time step 
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element matrices 

body force 

wave speed 

heat capacity 

arbitrary constants 

elastic constitutive tensor 

total energy density 

internal energy density function 

arbitrary function 

Heaviside step function 

Hilbert spaces 

determinant of Jacobian matrix for 
the transformation from inertial to 
reference frame coordinate system 
thermal conductivity (scalar and 
tensor) 

Boolean connectivity matrix 

outward normal unit vector to n 

shape function 

heat flux vector 

heat generated per unit mass 

time 

absolute temperature 

absolute temperature at zero strain 

displacement 

arbitrary function 

reference frame displacement 

relative displacement 

spacial coordinate 

thermal expansion coefficient 

boundary of open set 

frequency 

increment 

linear strain tensor components 
test function 
wave number 

Lame elastic constants 
natural coordinates 
density 

Cauchy stress tensor 
relaxation time 
shape function 


The use of heat and stress to alter materials and 
material surfaces has been known for over 5000 years. 
Although many of the compositions, heat treating, and 
joining techniques have been lost, it is evident that 
the objective was to combine the proper material com- 
position at a temperature and stress with heating and 
cooling rates to produce a more durable material. This 
is still the objective today; however, we are rapidly 
approaching some limitations in heat transfer rates 
and induced stresses. These processes are very 
important to our understanding of subsequent wear in 
gas path seals for turbine engines, as well as many 
other rapidly cooled surfaces. 

Conventional Cooling 

High heat transfer rates are experienced in boil- 
ing from a cylinder which passes through the stable 
regimes of film, transition, nucleation, and incipi- 
ence boiling prior to passing to the conduction and 
convection regimes. Nukiyama (ref.l) appears to ha/e 
first recorded this phenomenon, but it was undoubtedly 
well known to the heat-treating practitioners as a 
method of controlling material strength (eg. use of 
ice, salts, and oils). Thejjeak nucleate flux for 
water is about 15.8xl0 5 W/m* (5xl0 5 Btu/hr-ft 2 ) at 
DNB (departure from nucleate boiling) under steady- 
state conditions. However, the rate of cooling (heat- 
ing) associated with the Leidenfrost (incipience) 
temperature (which varies with bath subcooling, the 
material etc.) ranges from the metastable condition 
(ref. 2) to over 4.5xl0 3 K/s, as measured for bubble 
growth in saturated water (ref. 3). Giarrantano 
(ref. 4) investigated transient boiling in fluid 
helium for heating rates from steady state to the 
order of lxlO 4 K/s. The peak nucleate fluxes in- 
creased an order of magnitude over the steady-state 
values. 

If one increases the body force component (favor- 
ably) through flow augmentation devices (ref. 5) or 
by rotation (refs. 6,7) the heat flux can be increased 
by a factor of perhaps 20 (i.e., to 3xl0 7 W/m 2 
(lxlO 7 Btu/hr-ft 2 ) under steady-state conditions). 

The heat fluxes at the throat of a high-pressure 
rocket engine are 0.9x10° to 3x10° W/rrr (0.3x10° to 
1x10° Btu/hr-ft 2 ) and represent some of the highest 



operational steady fluxes. As to the temperature rise 
time of the materials during the starting and shutdown 
transients, which last for several milliseconds, the 
temperature changes are usually less than lxlO 4 K/s 
(nominally 3xl(r K/s). 

Rapid Quench 

In rapid quenching techniques and formation of 
metallic glasses, Davies (ref. 8) cites cooling rates 
which are very large. In the quenching literature the 
surface temperatures are used as boundary conditions 
(Dirichlet formulations) rather than the heat flux 
(Neumann formulations). Ice water, salt water, and 
oil quenching are commonly used in heat treating and 
forming of bulk materials. The rates are variable and 
depend on the mass, with fluxes and cooling rates both 
of the order of lxlO 6 K/s. The rates for producing 
metallic glasses all appear to be larger than 100 K/s 
to effect vitrification, with rates greater than 
lxlO 4 K/s generally considered practical.* 

Calculated values for the glass transition in metals 
lie between 3xl0 4 and lxlO 6 K/s (ref. 8). These 
rates are usually associated with at least one 
boundary in continuous motion, such as in turbine 
engine blade/seal rubbing. Laser glazing and splat 
cooling associated with plasma spraying are considered 
co produce the highest 

rates, 1x10 K/s and 1x10^ to 1x10^ K/s, respec- 
tively (ref. 8). 

For splat quenching of 40Fe-40Ni-14P-6B, with 
h = 2.5 MW/m — K and T^it — Tgiass =: 600 K 
(ref. 8), the average heat flux would be 1.5 GW/m 2 
(5x10° Btu/hr-ft 2 ). Thus although the cooling rates 
are very high, the heat flux requirements appear to be 
similar to those known to rocket designers. 

Mechanical Rubs 

In seal rubs, grinding, and wear phenomena, the 
problem is further complicated by the simultaneous 
application of a high normal load along with the 
moving boundary condition. Such problems may be 
classified as the Blok type (ref. 9): Blok studied 

the motion of a finite-width block sliding over a 
semi-infinite half space. Extension include dry rub 
geometries (turbine blades, seals, bearings, and 
brushes (refs. 10-12), and third-body effects (e.g., 
lubricant due to melt or conventional lubricants) are 
discussed by Godet et al. (ref. 13) and Braun et al. 
(ref. 14). 

The duration of the event is very short and more 
explosive in nature than either conventional cooling 
or rapid quenching; furthermore ablative losses can 
and do occur. Marscher (ref. 11) estimates the heat 
flux per unit of true contact area (including asperi- 
ties) to be 4xl0 3 Btu/in -s or 6.5 GW/m 2 (2xl0 g 
B tu/hr-f t** ) . Such high fluxes for even very short 
times are usually destructive because the stresses in 
the materials exceed the failure stress and surface 
cracking ensues. 

Objectives 

Although transient thermostress computations have 
been carried out for rocket engine channels and mech- 
anical rubs, the analyses are based on steady- and 
pseudo-steady-state postulates. In this paper, we 
develop an explicit transient formulation of the 
coupled thermal and mechanical equations, including 
wave effects, which can be used to predict stresses 
and thermal profiles under conditions of very rapid 
cooling or heating as associated with seal rubs, the 
formation of metallic glasses, other equivalent heat 
treat processes and rocket engine channels. For this 
information, boundary constraints of either the 
Lagrangian or Eulerian type can be handled. 

♦Exceptions occur e.g., metallic glass spherules of 
55-Au 22.5Pb-22.5Sb form at 10 2 to 10 3 K/s, ref. 32. 


MODELING AND MOTIVATION 

Consider, for example, problems with transient 
energy input and pressure loading of a fluid with 
deformable boundaries. The equations of motion for 
fluid dynamics (usually Eulerian in form) possess 
boundary conditions which are Lagrangian in form 
(usually associated with solid mechanics). The stan- 
dard methodology is to decouple the equations and 
introduce steady- or quasi-steady-state techniques. 

In this section we develop the equations and criteria 
for predicting the coupled thermomechanical behavior 
under very rapid and more conventional transients. 

The governing equation for the combined thermo- 
mechanical system is based on four conservation laws: 
conservation of mass, linear momentum, angular momen- 
tum, and energy. The form of the partial differential 
equations that represent these conservation laws'de- 
pends on the coordinate system and the kinematic 
description used. In this derivation the governing 
equations are constructed by using an arbitrary moving 
reference system as defined by its displacement U a , 
see figure 1. 

u = U + w 

~g ~ 

When Ug is zero, the description is the usual 
Eulerian formulation (stationary observer); when 
Da is the same as the material displacement u, the 
description reduces to the Lagrangian formulation 
fi«oving observer). 



FIGURE 1 . REFERENCE FRAME 


To complete the formulation, three sets of equa- 
tions are required in addition to the conservation 
principles: a constitutive set, an equation or set of 

equations of state, and a heat diffusion law or laws. 
However, to specify these equations requires a 
specific material and, since we wish to maintain gen- 
erality, these equations will remain arbitrary. 

The, use of this general description will allow the 
convenient numerical solution of problems involving 
large displacements in part of the mesh such as in 
sliding contacts (ref 15) and production of high- 
strength, corrosion-resistant ribbon materials (ref. 
16). Ribbons for transformer windings or pulverized 
for powder metallurgy and general ribbon materials of 
practical interest such as Fe, Ni, Co, Au, Al com- 
plexes are alloys which have predicted critical cool- 
ing rates for vitrification of the order of 3xl0 4 to 
1x10 K/s. The equations can be useful in problems 
involving large material flows such as in the analysis 
of continuous casting of steel (ref. 17). The formu- 
lation will also be useful in problems involving the 
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development of high-speed rubber products such as 
tires and crawler tracks where the speeds are lower 
but the applied loads are high. 


3 

3 1 



pu da 



pfc da 


GOVERNING EQUATIONS 
Conservation of Mass 

In this formulation the mass of any body is as- 
sumed to be fixed. The change of mass due "to chemical 
reactions is ignored. The change of mass in any open 
region a is then given by the mass flux into the 
region. Expressing this statement in integral form 
yields 



where p is the mass density, r is the boundary of 
the open region, and w is the velocity vector with 
Cartesian components wj of the mass with respect 
to the region a. 

The inertial velocity of the material jj can be 
decomposed into two components: the velocity of the 

material with respect to the reference frame, and the 
velocity of the reference frame w with respect to 
the inertial frame Ug (fig. 1). ~The vector iden- 
tity relating these three velocities is 


U = U g + w (2) 

In most previous work using an Arbitrary 
Lagrangian-Eulerian description, Hirt, Noh, Donea, and 
Belystchko (refs. 18-21) used the grid velocity as the 
independent variable describing the reference frame. 
When the grid velocity is used, the Eulerian descrip- 
tion can be obtained by setting the grid velocity to 
zero. The Lagrangian description, on~the other hand, 
can only be obtained by setting the grid velocity to 
the material velocity, which is unknown apriori. The 
preference of grid velocity over relative velocity is 
due to the application of the previous researcher's 
formulations to fluid problems where the Eulerian 
description is preferred. 

In the heat treatment problems the bulk material 
velocities are either known, as in the sliding-contact 
problem, or can be reasonably estimated, as in the 
continuous-casting process. The motion of the bound- 
ary of the material is not known apriori; therefore 
the description that most easily reduces to the 
Lagrangian description is used. Since equation (1) 
is true for an arbitrary boundary, by applying the 
divergence theorem and the Leibniz rule, the conserva- 
tion law can be written in differential form as 

~jr (pU) + 2 * (pfi)J = 0 (3) 

where J is the determinate of the Jacobian matrix. 

By using the identity due to Euler 



Ji - z • n) dr 


0 


( 6 ) 


The first term is the change of momentum within the 
volume a, the second term is the contribution from 
the body force vector j), the third term is the momen- 
tum flux out of the region, and the fourth term is the 
momentum on the body from the boundary traction. By 
converting equation (6) to a differential form in a 
similar manner as used in the mass conservation equa- 
tion, the following differential equation results: 


pu - pb - 2 • 2 + p(w 


S)i l = o 


(7) 


Conservation of Angular Momentum 

It can be easily shown that the conservation of 
angular momentum is satisfied only if the stress ten- 
sor is symmetric (which is the usual case). 

Conservation of Energy 

T H e conservat i o n of energy for a thermomechanical 
system with an arbitrary moving reference frame can be 
expressed in integral form as 


(2 • n) • & dr + 

/ pJj • fl dn - 


■'a(t) 

- / q • x) dr + 

J ro dQ = If 

■'aft) 3t 

J T - 

By splitting the energy density E 


'n(t) 


da 


( 8 ) 


into internal 

energy density & and kinetic parts and by using a 
similar procedure as in the derivation of the other 
conservation laws, a differential form of the conser- 
vation of energy can be written as 


a : £ - • if - s • q + pr = (9) 

To complete the governing equations, a constitu- 
tive equation, an equation of state, and a heat 
diffusion law must be specified. For example, 
restrictions of the form of these equations can be 
rationally shown as done in the works of Carlson, 
Trusedell, Noll, Gurtin, and Eringen (refs. 22-26). 
However, complete specification of these relations 
requires experimental information on a particular 
material. 

As an example, the constitutive relation for a 
linear thermoelastic solid is 



equation (3) can be rewritten in our final form as 

p + p 7 • J + vp • w = 0 (5) 

Conservation of Linear Momentum 

The conservation of linear momentum in an arbi- 
trary moving open region a with boundary r can be 
written in integral form as 


a . . — C. • ■ -j £ i •» 
ij i jk 1 kl 


*ij (T - V 


( 10 ) 


and the equation of state is given by 
e ij c ij + oCT + T 0°ij e ij 

When these relations are used with the classic Fourier 
law of heat diffusion 


1J 3 Xj 


( 12 ) 
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and isotropy is assumed, the energy equation win 
yield the linear coupled heat conduction equations as 
given in Boley (ref, 27). 

^M3x + 2uM 0 e kk =^(kfy (13) 

where x and u are the Lame elastic constants. 

However, in problems involving very rapid tempera- 
ture changes and large heat fluxes, which characterize 
the heat-treating processes, the assumption of linear 
behavior often leads to erroneous predictions. 

To develop useful numerical tool to investigate 
these problems, a formulation that is easily adapted 
to various nonlinear constitutive and state equations 
must be used. We therefore develop a finite element 
approximation to the field equations without assuming 
a form of the constitutive relations or a form of the 
equation of state. A general form of the heat diffu- 
sion equation is assumed as described in the following 
section. 

Heat Diffusion Equation 

The most valid first-order relation which exists 
between the heat flux vector and the temperature 
gradient is 


are assumed to be arbitrary, the explicit method pro- 
vides an effective formulation. 

The limitation for explicit methods is that the 
time increment used must be small enough to keep the 
integration stable. In rapid heat-treating problems 
this restriction is not excessive. The stability of 
our finite element formulation is investigated in the 
section Analysis of the Discrete Equations. 

By multiplying the differential forms of the con- 
servation equations by a test function n and inte- 
grating over the domain of the problem, the following 
weak forms for the conservation principles results: 

Conservation of Mass 


Find peH°, ui, wi e H 7 V n e such that 



Conservation of Linear Momentum 


q i 


-k 


il_ t £L 

3X i ij 3Xj 


(14) 


When this equation is used in a purely thermal problem 
with a fixed reference system, the energy equation be- 
comes the parabolic heat conduction equation 


pCT = kv 2 T (15) 

This equation, although adequate for slow processes, 
predicts an infinite propagation velocity for a finite 
heat flux. One method to resolve this error is the 
modification to equation (14) proposed by Maxwell and 
Grad (refs. 28, 29) 
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— - k — 
ax i ij axj 


T 
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3t 


(16) 


We use a modified form which is more convenient to 
incorporate into the finite element formulation: 



This diffusion law in a purely thermal problem will 
yield an energy equation that is of the form of the 
telegraph equation which was studied by Weymann, 
Baumeister, and Hamill (refs. 30, 31). 

Further heuristic generalization involving powers 
of the thermal gradients and multiple relaxation times 
can be made but are not considered in the proposed 
finite element formulation. 

FINITE ELEMENT FORMULATION 

Approximate solutions to the field equation will 
be developed by applying the finite element technique 
in space and a simple explicit integration method in 
time. The use of the explicit temporal integration 
method permits the transient solution to proceed with- 
out the need to solve simultaneous systems of non- 
linear equations. The behavior of current techniques 
for solving systems of nonlinear equations depends on 
the properties of the equations to be solved. Since 
the equations for the state and constitutive relations 


Find u-j, w*j e H^, p, a - jj e f n c such that 



Conservation of Energy 

Find & 9 T e H*, p, a-jj, e-jj e 


/.i 


n Vij ' Pn 3X 


_ .» _ in-!" i 

j j ,x jL‘ 


lP tn eH 1 such that 

lL + £ il_ 

3Xj K ij 3X. 


+ T ( k IxT + k ij 1 xt)] + n ° r “ pn ^ j da “ 0 (20) 

We will now develop a finite element approximation 
to the weak form of the conservation equations for 
two-dimensional problems. All functions will be 
approximated in space by assuming the form of the 
functions on each element to be given by 

jy(*,t) = N I (^)U J (t) (21) 

where the majuscule subscripts indicate the node num- 
ber and N are the shape functions defining the type 
of element. For a four-node bilinear quadri lateral 
element, shown in figure 2, the shape functions used 
to approximate H l functions in natural coordinates are 


N x = * (1 - C)(l - 4>) 
N 2 =| (1 + C)(l - 4-) 
N 3 = j (1 + S)(l + 4-) 

N 4 =| (1- 5)(1 + +) 


(22) 
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FIGURE 2. GENERIC ELEMENT 

The shape functions used to approximate H° functions 


are 


*!-[!- H(t)][l - HU)] 


» 2 = H ( C) Cl - HU)] 


(23) 


*3 = H( ?)H( h,) 

» 4 = [1 - H(«)]H(*) 

Using the shape functions in the weak form of the mass 
conservation equation results in tl^ following semi- 
discrete system of equations: 


B JI*J + A IJ P J " 0 


where 
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*K N P da L p j 
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(25) 
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The finite element approximation to the momentum equa- 
tion is obtained from the weak form in a similar 
manner. 

G ij u jk - Q ik + R u u Jk + s ik < 27 > 

Here the matrices G, Q, R, and S are defined by 
e 


3 IJ 


s * 4 


da L pj 
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(31) 


The finite element approximation to the energy equa- 
tion is 


G l/j = D I 

where 

e 


(32) 
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aN . 
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J 1 J 1 J 3 X, 


(33) 


All first-order time derivatives will be approximated 
by a forward difference operator 


;i f 1+1 - f 1 
At 


(34) 


and all second-order time derivatives will be approxi- 
mated by a central difference operator. 


f 1 - (Atf 2 (f 1 + 1 - 2f 1 + f 1 " 1 ) 


(35) 


When the temporal difference operators are applied 
simultaneously, the temporal integration scheme is 
always unstable. To mitigate the instability of the 
discrete system, a staggered approach to the temporal 
integration is used. The properties of the discrete 
equations are analyzed in the next section. 

ANALYSIS OF THE DISCRETE EQUATIONS 

The stability of the preceding discrete equations 
is investigated by examining several model equations. 

The first model equation is the Lagrangian formu- 
lation for one-dimensional elastic solids. The dif- 
ference equations for this system are 
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* r(i 1 * 2T i * - ii - <ii) 

* l '( 2 TJ n- T Ll- T ;u )- 0 

To study the stability of this system, it will be 
linearized by considering linear perturbations of 
u, and T. 


(38) 


p» 


TOTAL 

P 


(39) 


0 TOTAL 


0 o -0 


(40) 


One such approach is to advance one variable at a 
time. For example, advance the density by the standard 
forward difference operator; then, in the momentum 
equation, use the new values for the density to cal- 
culate the right side and solve for the linear dis- 
placement. The internal energy (temperature) can be 
advanced by using the new values for both the density 
and the displacements to compute the right-side vector. 

The linearized difference equation for this tem- 
poral integration method when applied to the linear 
problem can be shown to be stable when 



The second model problem is that of linear thermo- 
elasticity with a motion of the reference frame. By 
a similar analysis the conditions for a stable solu- 
tion are 


T T0TAL = T o + T (41) 

where the last term in each expression is assumed to 
be small. The linearized difference equations are then 

( p J+ i - p J i + 6p J+1 - 6p J + - p J . ) 

at\ M n-l p n-l .. M n M n M n+1 w n+l / 

+ 4 °o ( &n+l - “n-l) - 0 
( p n+l " p n-l ) = 0 


(42) 
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TF p o 


- 3c 


(43) 


« <(f y* 

(51) 

4t < :i r(r)( Ax)2 

(52) 

At < — 

(53) 


iWi 


The last restriction can be eliminated by intro- 
ducing upwinding. One way to accomplish this is to 
postmultiply all transport terms in the discrete equa- 
tions by an amplification matrix (ref. 22). 

SUMMARY 


It °o c [ T £i - t L * 4 ( T f 1 - T 0* T » : i - '«] 

* * - fj-i - ’in) * r T 0 ( a n - ‘jn - l ii) 

*r“o( T J*i* 2T J‘ T "-i)-° (M) 

Exact solutions to the linear difference equations in 
the form 

exp (iatyj + K nax) (45) 

ujj = C 2 exp (iatyj + K nax) (46) 

= C 3 exp (iatyj + <nax) (47) 

can be found; for the difference equations to be 
stable, the modulus of 

e 14tY , i = vCT (48) 

must be less than 1. Since this is not the case when 
the central-forward difference method is applied, a 
modified integration approach that is still explicit 
in form must be used. 


1. A formulation for transient heating which is 
independent of constitutive laws and the equation of 
state has been developed. The user selects a 
material, the appropriate constitutive law, and the 
equation of state; the numerical solution will then 
provide the thermomechanical effects for that material. 

2. The formulation applies to systems which re- 
spond to very rapid thermomechanical loadings such as 
occur in gas path seal rubs and in forming vitrified 
products where the quenching rates are greater than 
lxlO 4 K/s. It can also be applied to more 
conventional thermomechanical loadings but may not 
represent the most efficient technique. 

3. This formulation includes transport effects 
such as can occur in sliding contacts (e.g., bearings, 
seals, continuous casting, ribbon production, rubber 
products, and some types of heat treating). The cou- 
pled conservation laws predict two transport terms. 

The conventional heat transport iwi dT/dx and a 
strain energy transport term which is not included in 
the conventional governing equations for heat transfer 
(ref. 9). Studies need to be carried out to verify 
the magnitude of these terms. 

4. The stability analysis, although based on a 
linearized model, is effective because the estimate 
provides the maximum disturbance. Experience has 
shown systems with smaller perturbations to be 
stable. The time-step restrictions calculated from 
linear stabililty analysis can be used to estimate 
acceptable time steps in nonlinear problems. 

Experience indicates that from the maximum value 

of each coefficient (e.g., sound speed and thermal 
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diffusivity, etc.) the linear theory predicts reason- 
able values for the time increments in the nonlinear 
problem. These estimates can be used to judge the 
applicability of this formulation to a particular 
problem. 
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